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Abstract 


Global knowledge on the consequences of changes in vegetation and land 
use in the phases of the hydrological cycle and soil quality ¡is still incipient. 
This task ¡is particularly challenging ¡in tropical regions, where 
deforestation rates reach their highest reported levels and monitoring is 
limited. Such is the case of Mexico, where changes in coverage and water 
erosion affect more than half of the national territory. In this context, the 
SWAT model was used to study water erosion associated with variations 
in flow conditioned by the trajectory of change in vegetation cover in six 
periods, between 1986 and 2018, in two sub-basins of the Papaloapan 
River. Based on the combined evaluation of three efficiency and two 


uncertainty statistics, six of the ten calibrated periods showed good fitting 
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in discharge simulation with KGE > 0.70 and P-factor > 0.70. In the 
medium term, the results showed a decrease in sediment yield in the two 
study sub-basins, possibly related to the recovery of the forests in the 
basins' headwaters. In spite of the good results of the simulations, the 
detailed analysis of the sources of uncertainty in regions with limited 
monitoring is as important as the validation with good-quality measured 
data, in order to improve the performance and reliability of the models 


and to support adequate decision-making. 


Keywords: Erosion, flow, SWAT, deforestation, land use change, tropics, 


efficiency, uncertainty. 


Resumen 


El conocimiento global de las consecuencias del cambio de vegetación y 
uso de suelo en las fases del ciclo hidrológico y la calidad del suelo aún 
está en ciernes. Esta tarea es particularmente desafiante en las regiones 
tropicales donde las tasas de deforestación alcanzan los niveles máximos 
reportados y el monitoreo es limitado. Tal es el caso de México, donde el 
cambio de cobertura y la erosión hídrica afectan a más de la mitad del 
territorio nacional. En este contexto, se usó el modelo SWAT para estudiar 
la erosión hídrica asociada con las variaciones en el caudal condicionadas 
por la trayectoria de cambio de la cobertura vegetal en seis periodos entre 
1986 y 2018 en dos subcuencas del río Papaloapan. Basados en la 
evaluación combinada de tres estadísticos de eficiencia y dos de 


incertidumbre, 6 de los 10 periodos calibrados mostraron buenos ajustes 
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en la simulación de la descarga con KGE > 0.70 y P-factor > 0.70. Los 
resultados mostraron a mediano plazo una disminución en la producción 
de sedimentos en las dos subcuencas estudiadas, posiblemente 
relacionada con la recuperación de los bosques en las cabeceras de las 
cuencas. A pesar de los buenos resultados de las simulaciones, el análisis 
detallado de las fuentes de incertidumbre en regiones con monitoreo 
limitado es tan importante como la validación con datos medidos de buena 
calidad para mejorar el desempeño y la confianza de los modelos y 


soportar de forma adecuada la toma de decisiones. 


Palabras clave: erosión, caudal, SWAT, deforestación, cambio de uso de 


suelo, trópicos, eficiencia, incertidumbre. 
Received: 04/06/2020 
Accepted: 29/01/2022 


Introduction 


Water erosion is a naturally occurring and anthropogenic process of the 
earth's surface, whose magnitude on a global scale is estimated to be in 
the range of 20-30 gigatons per year (FAO € ITPS, 2015). The highest 
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rates of deforestation occur in the tropical regions of the Global South, 
(Borrelli et a/., 2017) promoted by the expansion of the agricultural 
frontier, the extraction of resources and/or growing urbanization Montes- 
León, Uribe-Alcántara, 8 García-Celis, 2011). During the last decade, 
transitions from forest cover to other land uses increased at a rate of 0.46 
% in Latin America and the Caribbean (FAO, 2011) and of 0.3 % in Mexico 
(Blackman, Goff, € Rivera-Planter, 2018); both rates exceed the world 
average of 0.13 % (FAO, 2011). Compounded by human activity, water 
erosion in the tropics is exacerbated due to high-intensity seasonal rainfall 
regimes (Rangel, Jorge, Guerra, € Fullen, 2019). Thus, more than 60 % 
of Mexican soil presents some level of erosion due to water degradation, 
which worsens at an average annual rate of 0.5 % (Borrelli et a/., 2017). 
These serious water erosion processes become especially acute in the 
stages subsequent to changes in cover and land use (Cardoza-Vázquez et 
al., 2007; Bolaños-González et a/., 2016). 


Changes in coverage affect environmental services, the water cycle 
and the soil at multiple spatial and temporal scales, depending on the 
physiography of the basin, climate variability, and the trajectory of use of 
the studied territory (Martinez, Guyot, Filizola, 8 Sondag, 2009; Arheimer 
8: Lindstróm, 2019). For example, vegetation and land use are important 
factors in the hydrological behavior of basins, since they determine 
surface and subsurface runoff generation processes  Arheimer 8 
Lindstróm, 2019; Couto et a/., 2019; Liu, Zhang, Kou, € Zhou, 2017). At 


a regional scale, the volume of runoff and transported sediments 
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increases, causing morphological alterations to the river network and the 
hydrological regime of the basins, and there is a higher risk of flooding in 
the plains (Chotpantarat 8: Boonkaewwan, 2018; Wang, Ge, Yu, Wang, € 
Xu, 2017; Yang et al., 2018). 


In tropical countries, hydrological monitoring is performed with 
scanty measurement networks (Vórósmarty, Lévéque, 8 Revenga, 2005), 
resulting in limited historical records (Carvalho, 2019; Delmas, Cerdan, 
Cheviron, 8: Mouchel, 2011; Warrick, 2014). The convergence of historical 
data scarcity, the need of information for decision-making, advances in 
computing technologies, and the production of scientific knowledge 
contribute to the application of hydrological models as the most widely 
used tools for research, planning, and environmental management 
(Fakhri, Dkohaki, Eslamian, Fazeli-Farzani, 8 Reza-Farzaneh, 2014; 
Jajarmizadeh, Harun, € Salarpour, 2012; Joórgensen 8 Fath, 2011). 
Similarly, geographic information systems have enhanced the 
management of large amounts of information, yielding advantages such 
as the display of spatially distributed results and the generation of 
scenarios (Alatorre € Beguería, 2009; Goodchild, 1996). 


The Soil and Water Assessment Tool (SWAT) is a physically based 
conceptual model that estimates water, sediment, nutrients, pollutants, 
and pathogens at basin level from spatially distributed data on 
topography, soil, vegetation and land use, management practices, and 
climate (Douglas-Mankin, Srinivasan, € Arnold, 2010). The operation of 


this model is defined by the principles of mass and momentum: the former 
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is related to the balance of water, sediments, nutrients or contaminants, 
and the latter focuses mainly on the displacement of these flows through 
the drainage network of the basin. The main areas of application of the 
SWAT model include water quality in relation to non-point source 
contamination, sediment production, runoff simulation, and the 


hydrological impacts of environmental change (Wang et a/., 2019). 


In the SWAT model, the dynamics of the variables of interest is 
modeled from spatially distributed physiographic, climatic, and anthropic 
input data, so that it is useful in ungauged basins or in the evaluation of 
long-term scenarios based on variations of some of the input variables 
(Neitsch, Arnold, Kiniry, € Williams, 2011). In this case study, it will be 
based on changes in coverage and land use. The use of the SWAT model 
in countries with limited monitoring is widespread. Across regions of East 
Africa, Guzha, Rufino, Okoth, Jacobs and Nóbrega (2018) identified 15 
studies that used SWAT to assess the impact of land use on surface runoff, 
discharge, and baseflow. For their part, Tan, Gassman, Srinivasan, Arnold 
and Yang (2019) reviewed 126 cases that studied different issues of land 
use change, climate change, management practices, and/or water quality 
with SWAT in Southeast Asia. 


In the Papaloapan River basin, Bello, Gómez, Magaña, Graizbord 
and Rodríguez (2009), and Ruiz-Fernández et a/. (2014) pointed out that 
siltation problems in the system of lagoons and wetlands of the lower 
basin are associated with a large amount of sediment transported by 


tributaries upstream during the last forty years. Pérez-Vega and Ortiz- 
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Pérez (2002) established a positive relationship between the increase in 
the runoff of the sub-basins of the San Juan and Tesechoacán rivers and 


the increase in deforestation in the Papaloapan River basin between 1973 
and 1993. 


From the standpoint that environmental processes are historical 
—as a product of a historical legacy but also the result of a combination 
of variables in the present Reid 1998; Verstraeten Lang, € Houben, 
2009)— and susceptible to research with available resources, the SWAT 
model was used to reconstruct the trajectory of the potential erosion 


associated with land use change with the following specific objectives: 


1. To estimate water erosion and flow variations with respect to changes 


in land use and cover from 1985 to 2018 in two sub-basins of the 
Papaloapan River. 


2. To estimate the spatial distribution of the erosion potential in the study 
area. 


Study area and databases 


The sub-basins of the San Juan (SJN) and Tesechoacán (TES) rivers 
constitute 26 % of the 46 720 km? of the total area of the Papaloapan 
River basin, which is located on the slope of the Gulf of Mexico between 
19% N 972 W and 17% N 942 w (Figure 1). 
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Figure 1. Location of the study sub-basins (upper right) with elevation 
(left) and data grid from the Global Land Data Assimilation System 
(GLDAS) (lower right). 
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About 55 % of the study area has slopes that are less than 15 % 
steep and 45 % of the remaining area has slopes steeper than 15 % (FAO, 
2009). Using the Kóppen classification modified for Mexico by García 
(2004), three climate groups are identified in the region: (A) tropical 
rainy, in the plains close to sea level; (B) dry with summer rains and (C) 
humid temperate in the high basins (up to 3 400 masl). Twelve types of 
vegetation and land use were identified, grouped into five classes: 
Temperate forests, tropical forests, scrub and chaparral areas, pastures, 
and agriculture (Alavez-Vargas, Birkel, Corona, 8 Breña-Naranjo, 2021) 
(Figure 2). 


Tecnología y ciencias del agua, ISSN 2007-2422, 


Open Access bajo la licencia CC BY-NC-SA 4.0 14(4), 349-414, DOI: 10.24850/j-tyca-14-04-08 
(https: //creativecommons.org/licenses/by-nc-sa/4.0/) 


"o Check for ul dates 
open (o) access ¡0 E ! 


Tecnología y 


CienciaszAgua 


100% 
90% 90% ” 
80% 80% 
70% 70% 
60% 60% 
50% 50% 
40% 40% 
30% 30% 
20% 20% 
10% 10% 

0% 0% 


100% 


1986 1993 1998 2003 2011 2018 1986 1993 1998 2003 2011 2018 
Land cover 
Temperate forests Tropical forest Shrublands and scrublands Grasslands Agriculture | Sub-basins 
MM Oak forest MA Tropical rainforest HN] Shrubland Cultivated grassland PB Irrigated crops 
[A Coniferous forest Tropical dry forest PY Crasicaule scrubland  Induced grassland PB Annual crops 
HA Cloud forest ES Rosetophyllous scrubland Other aspects 


—— Drainage network 


Figure 2. Vegetation and land use from 1986 to 2018, geographical 


representation and percentages by study sub-basin. 


The Air Shuttle Radar Topographic Mission (SRTM) digital elevation 
model, with a spatial resolution of 90 m, was used in this study (Jarvis, 


Reuter, Nelson, 8 Guevara, 2008). This raster information was reclassified 
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into slope ranges, based on the geomorphology of the basin and FAO's 
classification (FAO, 2009) in order to create two classes: 0-15 % and 
steeper than 15 %. The vector sets corresponding to the drainage network 
(INEGI, 2010) come from the thematic geospatial information produced 
by the National Institute of Statistics, Geography and Informatics (INEGI) 
of Mexico. The edaphological information required by the model was 
derived from the vector sets related to the edaphological units (INEGI, 
2007) and the edaphological profiles (INEGI, 2013), in addition to two 
software packages for estimating other characteristics not available in the 
sources previously consulted: using the Soil Water Characteristics 
software (Saxton 8 Rawls, 2006), the saturated hydraulic conductivity, 
the apparent density, and basic infiltration were assessed, whereas the 
hydrological soil groups were identified using the NumCur 10 IE software. 
The types of vegetation and land use were classified according to the 
procedure described in Alavez-Vargas et al. (2021) from Landsat images 
downloaded from Google Earth Pro (Wuthrich, 2006), processed in QGIS 
vV. 3.4.5 Madeira (QGIS Development Team, 2019), for the years 1986, 
1993, 1998, 2003, 2011, and 2018 (Figure 2) and subsequently 
reclassified according to the nomenclature used by the model. Daily (Li, 
Beaudoing, € Rodell, 2018) and tri-hourly (Beaudoing 8 Rodell, 2015; 
Beaudoing €: Rodell, 2016) temperature, precipitation, wind speed, and 
solar radiation data, with an approximate spatial resolution of 25 km?, for 
the period 1980-2018, were obtained from the Global Terrestrial Data 
Assimilation System (GLDAS) portal (Rodell et a/., 2004) in a 63-point 


grid over the study area (Figure 1). Finally, the daily flow measurements 
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from 1985 to 2014 were taken from the National Surface Water Data Bank 
(BANDAS for its acronym in Spanish) (Conagua, 2017). The products, 
resolutions, and sources consulted to obtain the data of the ten 


biophysical and climate variables collected are summarized in Table 1. 


Table 1. Characterization of input data. 


Product and temporal 


Variable Resolution Source 
resolution 
Digital 
Elevation SRTM, raster format 90m Jarvis et al. (2008) 
Model 
Water Hydrographic network, vector Scale 
INEGI (2010) 
network format 1:50,000 
Scale 
Edaphological units INEGI (2007) 
1:250 000 
Edaphology 
Scale 
Edaphological profiles INEGI (2013) 
1:250 000 


Landsat images, raster format, 
Vegetation years 1986, 1993, 1998, 2003, 30 m Wuthrich (2006) 
2011, 2018 
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1) Check for updates 


Variable 


Product and temporal 


resolution 


Resolution 


Source 


Temperature 


GLDAS Noah Land Surface Model 
14 V 2.0 
Period: 1980-01-01 a 2010-12-31 


GLDAS Noah Land Surface Model 
14 V 2.1 


Three hours, 
0.259 x 0.250 


Three hours, 


Beaudoing y Rodell 


(2015) 


Beaudoing y Rodell 


0.25% x.0,259 (2016) 
Period: 2011-01-01 a 2018-12-31 
Precipitation | GLDAS Catchment surface Model 
A Daily Li et al. (2018) 
Wind 0.250 x 0.250 
velocity Period: 1980-01-01 a 2010-12-31 
Solar GLDAS Noah Land Surface Model 
vada 3 horas Beaudoing y Rodell 
radiation 0.250 x 0.250 (2016) 
Period: 2010-01-01 a 2018-12-31 
Location of the hydrometric 
Hydrometric stations Garro (Tesechoacán 
] ] Conagua (2017) 
stations River) and Cuatotolapan (San Juan 
River) 
Measured BANDAS daily records ] 
Daily Conagua (2017) 
flow Period 1985-01-01 a 2014-12-31 
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Methods 


SWAT model configuration 


SWAT modeling consists of five stages (Figure 3). Initially, the digital 
elevation model and the river network guide the delimitation of the study 
sub-basins. Subsequently, the variables with which the model defines the 
hydrological response units (HRUs) are incorporated: slope, soil type, 
vegetation type, and land use. These first two stages, therefore, focus on 
the physical delimitation of the units (sub-basins) and sub-units of 
analysis (HRUs) by which, in the next stage, the information of the input 
variables is grouped and the responses of the components of the water 
balance (Equation 1) and the variables of interest are calculated (erosion 
in this case, Equation 2) (Neitsch et a/., 2011). 
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Figure 3. SWAT modeling workflow: input data preparation and model 


Parameterization 


Scenario e 


application and calibration. 
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In the third stage, the model was fed with the following climate 
variables: precipitation, maximum and minimum temperatures, solar 
radiation, and wind velocity. Together, these variables represent the 
amounts of water and energy that govern the water balance. Specifically, 
these data are used by the model to solve the water balance equation 
(Equation 1) at HRU scale per sub-basin (Neitsch et a/., 2011): 


SW, = SWo + YXiz1 E > Qsurp 7 Ea = Wscep — Ou) (1) 


Where SW:: final soil water content (mm), SWo: initial soil water 
content on day ¡ (mm); t: time in days; Raay: precipitation on day ¡ (mm); 
Qsurf: surface runoff on day ¡ (mm); E.: evapotranspiration on day ¡ (mm); 


Wseep: percolation on day í¡ (mm), and Qyw: return flow on day ¡ (mm). 


In turn, erosion in HRUs by sub-basin is calculated with the Modified 
Universal Soil Loss Equation (MUSLE, Equation (2)) (Williams, 2012): 


0.56 
sed = 11.8 (Os X peak * ACA pru) X Kysie X Cusie X Pusie X LSysig X CFRG 


(2) 


Where sed: daily produced sediment (ton); Qsurf: daily surface 


runoff sheet (mm ha1); qpeax: peak flow (m3 s"1); areanru: HRU area (ha); 
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Kus.e: soil erodibility factor; Cus.e: plant cover and management factor; 
Puse: support practices factor; LSus.s: topographic factor, and CFRG: soil 


fragmentation factor. 


The result of the third stage is the initial simulation of the scenario 
in turn. This simulation is calibrated with data of the analyzed variables, 


as described in the following section. 


Initial simulation and calibration 


This study was carried out in QGIS, version 2.6.1 Brighton (QGIS 
Development Team, 2014) with the QSWAT extension, version 1.8 (Dile, 
Srinivasan, € George, 2019). The input data was prepared in the format 
required by the software for the calculation of the different variables 
involved in equations 1 and 2. The equations given by default by the 
model were used, (see Arnold et a/., 2011) except for evapotranspiration, 
which was calculated with the Hargreaves, Hargreaves and Riley (1985) 
method. The values of the physical parameters (¡.e., the elevation, the 
drainage network, the slope and the edaphological characteristics) were 
constant for the six periods assessed in each sub-basin. Vegetation covers 


and land use were specific for each scenario. The climate information 
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introduced had a daily resolution. Due to the absence of databases on soil 
conservation practices in the study area, the Pus. support practices factor 


was considered equal to 1, ¡.e., without conservation practices. 


For the period from 1986 to 2018, six models were carried out per 
sub-basin, corresponding to the periods defined by the vegetation and 
land-use scenarios, as shown in Table 2. The irregularity of the assessed 
intervals is due to the methodological decisions made to ensure the 
mosaic of images of the best possible quality to characterize the cover 


and land use in each scenario. 


Table 2. Periods of warming, modeling and calibration in the three 


studied sub-basins. 


Vegetation and Monthly flow 
land use scenario dd dai calibration 
1986 1982-1986 | 1982-1983 1985-1986 
1993 1985-1993 | 1985-1986 1987-1993 
1998 1992-1998 | 1992-1993 1994-1998 
2003 1997-2003 | 1997-1998 1999-2003 
2011 2002-2011 2002-2003 2004-2011 
2018 2010-2018 | 2010-2011 2012-2014 
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For each modeling, a warm-up period of two years was defined, as 
suggested by Arnold et al. (2012). The modeling had a monthly time 


scale. 


Flow calibration was performed with 30 years of observed data 
(1985-2014) at the hydrometric stations of the San Juan and 
Tesechoacán sub-basins. The simulated erosion was not calibrated due to 


the absence of erosion data in the three study sub-basins. 


The sensitivity, uncertainty, and calibration analyzes of the models 
were carried out following the recommendations of Abbaspour, Vaghefi 
and Srinivasan (2017) with the SUFI-2 algorithm available in the SWAT- 
CUP software (Abbaspour, 2015). The SUFI-2 algorithm focuses on the 
search for the group of parameters with the best fit and least uncertainty 
regarding the data observed by means of sampling by the Latin hypercube 
method (Beven, 2012), with the aim of fulfilling two conditions: that more 
than 90 % of the data be located in the 95 percent prediction uncertainty 
(95PPU) area and that the proportion in the average difference between 
the limit values of the 95PPU and the standard deviation be less than 1 
(Abbaspour, Johnson, € van Genuchten, 2004). 


For each scenario and sub-basin, at least three iterations of 600- 
1000 simulations were carried out. The Kling Gupta Efficiency (KGE) 
criterion was used as the objective function, and the Nash-Sutcliffe 
efficiency statistics (NSE), the percentage of bias (PBIAS) and the 


coefficient of determination (R?) were also performed. 
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The Kling Gupta criterion evaluates the goodness of fit of the model 
by calculating the Euclidean distances between the correlation, the bias 
and the variability of the simulated value with respect to the observed 
value (Equation (3) and Equation (4)) as an alternative to the Nash 
Sutcliffe criterion, which takes as reference the difference of 1 minus the 
ratio of the mean squared error and the variance (Equation 5) to show 
the proportion of the variability of the observations explained by the 


simulation (Gupta, Kling, Yilmaz, 8, Martínez, 2009): 


KGE =1-ED (3) 


ED=J(r=D?+(a- 1)? + (8 - 1)? (4) 


Where ED is the Euclidean distance from the ideal point, r is the linear 
correlation coefficient between simulated and observed data, ais the 
variability given by the ratio between the standard deviations of the simulated 
and observed data, and fis the bias given by the ratio between the 


simulated and observed mean flows: 


n ES 2 
NSE = 1 Hlrsito) q M5E dE 


Z 
E ,(%o,+-Ho) 9 
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Where n is the total number of time steps, Xs,: is the simulated value 
at time t and uo and oo are the mean and standard deviation of the 


observed values. 


The bias percentage assesses the average tendency of the 
simulated data to overestimate or underestimate the observed data. The 
P and R factors are uncertainty statistics about the 95PPU prediction 
interval, where the P-factor refers to the percentage of observed data 
falling within the 95PPU prediction band, and the R-factor indicates the 
width of the 95PPU prediction band (Abbaspour, 2015), which the 


narrower, the less uncertainty it indicates (Abbaspour, 2005). 


From the initial best simulation and the sensitivity analyses, the 
software suggests a series of meaningful parameters to perform the 
calibration (Table 3). In each of the study sub-basins, flow calibration was 
carried out for each of the six scenarios, with a variable number between 
four and eight parameters (Table 4 and Table 5), depending on the 


sensitivity of each simulation to variations in the parameters involved. 
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Table 3. Parameters used in the discharge calibration of the three sub- 


basins. 
o Hydrologycal 
Acronym Description of the parameter 
component 
SCS runoff curve number for 
R_CN2.mgt Surface runoff 


V_ALPHA_BF.gw 


V_GW_DELAY.gw 


moisture condition II 


Groundwater delay (days) 


C2 | V_ALPHA_BF.gw | Base flow alpha factor (days 1) |] 


Groundwater 
Threshold water level in shallow 
V_GWQMN.gw 
aquifer for base flow (mm H20) 
Soil available water storage , 
R_SOL_AWC(..).sol Soil water 


capacity (mm H20/mm soil) 


R_ESCO.hru 


Soil evaporation compensation 
factor at HRU level 


R_ESCO.bsn 


Soil evaporation compensation 


factor at basin level 


V_GW_REVAP.gw 


V_REVAPMN.gw 


Groundwater revaporization 


coefficient 


Threshold depth of water in 
shallow aquifer for revaporization 


to occur (mm H20) 


Potential and real 


evapotranspiration 
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Table 4. Initial and final ranges of the parameters used in the flow 


calibration and efficiency statistics of the TES basin for each of the 


Acronym 


R_CN2.mgt 
V_ALPHA_BF.gw 
V_GW_DELAY.gw 
V_GWQMN.gw 
R_SOL_AWC(..).sol 
R_ESCO.hru 


R_ESCO.bsn 


1985-1986 
-0.1/ 0.11 
0.23/ 0.36 
98.03/ 192.86 
1.95/ 2.42 
-0.34/ 0.72 


-0.12/ 0.51 


modeled periods. 


1987-1993 
-0.16/ -0.12 

0.26/ 0.35 
58.75/ 95.52 
111.32/ 178.11 


0/ 0.06 


0.31/ 0.72 


Initial range (min/max) 


Period 

1994-1998 1999-2003 
-0.2/ -0.07 -0.2/ -0.12 
-0.01/ 0.95 -0.01/ 1.02 
-3.08/ 262.89 — 277.78/ 541.58 
-1.38/ 0.34 -9.17/ 144.7 
-0.34/ -0.14 0/ 0.77 
-0.41/ -0.14 0.27/ 1.61 


2004-2011 
-0.2/ -0.15 
-0.05/ 0.2 
20.21/ 152.76 
1.15/ 2.11 


0.02/ 0.48 


0.63/ 0.94 


See Table 3 for the description of the parameter represented by the acronym. 
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2012-2018 
-0.2/ -0.18 
0.29/ 0.4 
123.94/ 203.7 
0.53/ 1.2 


-0.88/ -0.41 


-0.02/ 0.26 
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Table 5. Initial and final ranges of the parameters used in the flow 


calibration and efficiency statistics of the SJN basin for each of the 


modeled periods. 


Initial range (min/max) 


No. Acronym Period 
1985-1986 1987-1993 1994-1998 1999-2003 2004-2011 2012-2018 
1 R_CN2.mgt -0.09/ 0.15 -0.09/ 0 -0.2/ 0.05 -0.2/ 0.06 -0.19/ -0.16 -0.2/ 0.04 
2 V_ALPHA_BF.gw 0.32/ 1.03 0.02/ 0.39 -0.11/ 0.66 -0.4/ 0.55 0.39/ 0.61 0.45/ 1.4 
3 V_GW_DELAY.gw 207.91/ 587.09 | -8.05/ 122.06 | -139.63/ 262.63 49.77/ 325.23 290.53/ 355.01 | -107.24/ 272.24 
4 V_GWQMN.gw 0.1/ 1.4 1.48/ 2.15 0.94/ 2.96 44,59/ 140.41 121.19/ 156.32 -0.65/ 1.15 
5 R_SOL_AWC(..).sol -0.57/ 0.33 0.51/ 0.66 0.08/ 0.71 
6 R_ESCO.bsn -1.3/ -0.45 0.45/ 0.67 0.25/ 0.8 
7 V_GW_REVAP.gw 0.16/ 0.21 
8 V_REVAPMN.gw 17.94/ 159.06 


See Table 3 for the description of the parameter represented by the acronym. 
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Results 


Land use change 


The area of the TES sub-basin covered by temperate forests increased 
from 16 % in 1986 to 26 % in 2003 and then declined again in 2018 
(Figure 2, graph of TES areas). For their part, the rain forests showed a 
50 % loss in their extension from 1986 to 2018, with a marked decline in 
2003 and, later, a recovery in the following 15 years. The areas occupied 
by chaparrals and rosetophyllous scrubs were less than 0.25 % in all 


scenarios. 


Variations in natural plant cover are directly related to variations in 
production practices. The recovery of temperate forests was mainly due 
to the discontinuation of rainfed agriculture in the highlands between 
1986 and 2003, while the decline in the period from 2003 to 2018 
indicates human practices that promote deforestation. In turn, the 
decreased extension of the tropical forest in the lower parts occurred due 
to the clearing carried out for the expansion of pasture areas (for livestock 
purposes), which came to occupy half of the sub-basin in 2003. In 
subsequent years, the abandonment of the grasslands reduced their area 


to occupy only 30 % of the basin, allowing the recovery of secondary 
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tropical forest vegetation in those places. From 1986 to 2003, the 
agricultural areas gradually reduced from 12 % to less than 1 % of the 


area. Later, in 2018, they expanded to occupy 20 % of the sub-basin. 


The San Juan sub-basin is the southernmost one and has the lowest 
altitude (< 1 000 masl), so there are less areas covered with temperate 
forests (<5 %, except in 2003 and 2011) than in TES (Figure 2, graph of 
SJN areas). On the contrary, tropical forests covered 55 % of the area in 
1993, followed by a degradation process that reduced it to a quarter in 
2003 and a recovery phase in the last two decades. In all scenarios, 


chaparrals represented less than 0.05 % of the area (>400 ha). 


Similarly to TES, the area lost to forests and tropical forests from 
1986 to 2003 was gained by grasslands, which by 2003 had covered 68 
% of the sub-basin and was later reduced to 39 %. As for agricultural 
areas, they have increased the occupied area from 10 % in 1986 to 25 % 
in 2018, with the lowest point reached in 2003. However, as discussed 
below in relation to the quality of the input data (see section The sources 
of uncertainty and the scope of modeling), the drastic changes registered 
in 2003 with respect to the previous (1998) and subsequent (2011) 
scenarios suggest the need to verify the results obtained for that year 


with a different set of images. 


Tecnología y ciencias del agua, ISSN 2007-2422, 


Open Access bajo la licencia CC BY-NC-SA 4.0 14(4), 349-414, DOI: 10.24850/j-tyca-14-04-08 
(https: //creativecommons.org/licenses/by-nc-sa/4.0/) 


o) 1) Check for updates 
OPEN ACCESS 
Tecnología y 


CienciaszAgua 


Hydrological modeling 


The results presented in this section and the next are the product of the 
parameterization of the models with the best fit obtained according to the 
efficiency statistics. Table 4 (TES) and Table 5 (SIN) show the initial 
uncertainty ranges given by default by the model in reference, to which 


the final values of each parameter for each sub-basin were adjusted. 


The graphs with the parameterized flow data with the best fit and 
the prediction intervals for each period are presented in Figure 4 and 


Figure 5, together with the respective efficiency and uncertainty statistics. 
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Figure 1. Flow modeling for 1986-2014 (calibrated) and 2015-2018 


(simulated) at the “Garro” hydrometric station of the Tesechoacán (TES) 


sub-basin. Each graph (section) corresponds to a modeled period (see 


Table 2). 
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Figure 2. Flow modeling for 1986-2014 (calibrated) and 2015-2018 


(simulated) at the "Cuatotolapan" hydrometric station of the San Juan 


(SIN) sub-basin. Each graph (section) corresponds to a modeled period 
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Regardless of the changes in vegetation and land use between the 
different scenarios by sub-basin, the same six parameters presented the 
best fit in TES (Table 4). 


The bias percentage indicates that the results are good when they 
are in the range of +15 % according to Moriasi et a/. (2007). Abbaspour 
(2015) suggests that, in discharge calibration, P-factor values >70 % are 
considered optimal. The same author points out that the R-factor value 
must be less than 1.5, and for this sub-basin all the periods presented 
values between 0.31 and 0.58. 


TES presents efficiency values in all the calibrated periods between 
0.64 and 0.90, except in the 1999-2003 period (KGE = 0.58, NSE = 0.14, 
R? = 0.35) (Figure 4). The poor performance of the 1999-2003 period is 
verified in the low P-factor value = 0.28, the comparatively high bias 
coefficient (PBIAS = 7.5), as well as visually (Figure 4d). The uncertainty 
statistics confirm the satisfactory adjustment in the 1994-1998, 2004- 
2011 and 2012-2018 periods (values of P-factor greater than 70 % in the 
three cases) and indicate that the 1985-1986 period is slightly below the 
performance threshold (P-factor = 0.63). The first and third periods are 
slightly overestimated, while the rest of the periods present a slight 
underestimation. Although the period of the simulated data from 1987 to 
1993 presented the best values of KGE, NSE and R? (> 0.80), only 40 % 
of the data observed in 1987-1993 coincide with the 95PPU interval 
(Figure 4b). 
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In SJN, surface runoff and groundwater, with one (R_CN2.magt) and 
three (V_ALPHA_BF.gw, V_GW_DELAY.gw, V_GWQMN.gw) parameters, 
respectively, explained the behavior of water yield from 1986 to 1998, 
while from 2003 to 2018, soil water (R_SOL_AWC(..).sol) and potential 


and real evapotranspiration had greater relevance (Table 5). 


Similarly to TES, in SJN (Figure 5), the 1999-2003 period shows 
the lowest efficiency values (< 0.54). In contrast, periods 1987-1993, 
1994-1998, and 2012-2018 showed the best performance values in the 
KGE (0.90, 0.78, and 0.89, respectively), NSE (0.81, 0.57, and 0.83, 
respectively), and R? (0.83, 0.62, and 0.85, respectively) statistics, as 
well as in the P-factor (0.82, 0.72, and 0.75, respectively) and R-factor 
(0.72, 0.72, and 0.75, respectively). Although the 1985-1986 and 2004-— 
2011 periods show good adjustments, the P-factor indicates a higher level 
of error with values of 0.54 and 0.57, respectively. 


Modeling potential erosion 


The potential erosion values estimated in this section correspond to the 
multi-annual averages of the six modeled flow periods. The average water 
erosion per scenario per sub-basin oscillates in the range of 69 to 158 ton 


ha"? year? in the Tesechoacán sub-basin and between 57 y 162 ton ha"? 
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year? in San Juan (Table 6). According to these figures, the volume 
eroded in Tesechoacán has decreased to 44 % of the volume estimated 


in 1986, and in San Juan to 38 % of the initial estimate. 


Table 6. Multi-annual average water erosion (ton ha"? year1) per period 


per sub-basin. 
Multi-annual average water erosion (ton ha"! year!) 
Sub-basin Period 
1985-1986 1987-1993 1994-1998 1999-2003 2004-2011 2012-2018 
Tesechoacán 158.4 73.83 108.36 81.1 80.46 69.17 


San Juan 162.28 67.17 73.86 100.42 56.77 61.98 


Figure 6 shows the spatial distribution of the erosive potential of 
HRUs by sub-basin, classified into four erosion categories: none to light, 
moderate, high, and very high (FAO, 1980). In TES, the greatest erosion, 
especially notable in 1986, occurred in the upper part of the sub-basin, 
where temperate forests are concentrated on steep slopes, followed by 
the middle basin, with moderate to high erosion, while the lower part 
reported light to moderate volumes in all scenarios. In the SJN sub-basin, 
on the other hand, the areas of greatest erosion are also located in the 
areas of highest elevation, although, during the different periods, small 


areas with very high erosion are observed in the middle and low areas. In 
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this sub-basin, a reduction in the area affected by more severe levels of 


water erosion is observed. 
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Figure 3. Distribution of the erosive potential by sub-basin according to 
FAO's classification (FAO, 1980). 
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In percentage terms (Figure 6, graph of TES areas), more than 50 
% of the TES basin shows moderate erosion during all the scenarios, 
except in 1986 and 1998, when the area with high erosion increased. 
Areas with light and very high erosion represented each less than 20 % 
of each scenario, on average. In the SJN sub-basin (Figure 6, graph of 
SJN areas), between 40 % and 55 % of the area reported erosion levels 
between 10-50 ton ha"! year, except in 1986; all other levels remained 


mostly below 20 %. 


Discussion 


Precipitation variability and observed flows 


For the 1986-2018 period, the average rainfall was 1,795 mm in TES and 
2,250 mm in SJN. Both sub-basins registered the lowest levels of 
precipitation in the 1994-1998 period, with up to 450 mm less (in SJN) 
than the above-mentioned average, coinciding with one of the periods 
(1997-1998), in which the El Niño-Southern Oscillation (ENSO) 


meteorological phenomenon had a greater impact on the national territory 
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(Escobar-Briones, Bonilla, Badán, Caballero, € Winckell, 2001). This 
decrease was observed in reductions of between 5 and 10 % of the 
average flow for the 1994-1998 period in the two sub-basins, with respect 
to the other periods. In contrast, precipitation for the 2012-2018 period 
in SJN exceeded more than 50 % the average of the 34 years analyzed, 
with atypical monthly maximum values above 500 mm in the rainy 
seasons of 2017 and 2018 (Conagua, 2018). 


The average flow observed for the six study periods ranged between 
151 and 196 m* s* in TES and between 208 and 251 m* s* in SIN. This 
variability in each basin is explained on two time scales. First, the 
seasonality of the annual precipitation is reflected in the amplitude of the 
monthly ranges of the values of flows observed and simulated for each 
period evaluated by sub-basin (Figure 4 and Figure 5). Second, the 
influence of ENSO is responsible for reductions of about 10 % in summer 
rainfall and increases of about 20 % in winter rainfall in the years 1986- 
1987, 1991-1992, 1994, 1997-1998, and 2002 (Magaña, Pérez, € 
Conde, 1998; Magaña, Vázquez, Pérez, € Pérez, 2003). 
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Land use change and sediment production 


From 1986 to 1993, TES and SIN showed a reduction of more than 50 % 
in the volume of erosion produced (Table 6), attributable to the increase 


in forest and rainforest cover (Figure 2). 


In the 1998 period, an increase was observed in the area with the 
greatest severity of erosion —from moderate (10 ton ha"! year?) to high 
(10-50 ton ha? year*)— in the middle zone of the TES sub-basin (Figure 
6). This exacerbation is mainly related to the expansion of grasslands, 
which, although they continued to expand until 2003 (Figure 2), no longer 
reflected the increase in average erosion of that period. In contrast, in 
the 2003 scenario, a decrease in the area affected by high erosion and a 
25 % lower sediment production is observed in TES compared to the 
previous period (Table 6). This decrease is possibly associated with the 
reduction of agricultural land use in favor of forest recovery, particularly 
in the upper basin, similar to what was found by Muñoz-Villers and López- 
Blanco (2008) in La Antigua, a basin of the same hydrological region, in 
the same period. However, the 2003 vegetation and land use scenario 
presented the most drastic changes in primary vegetation in favor of 
grasslands (Figure 2). In SJN, the expansion of grasslands is also marked, 


especially in the tropical forests, and this change translated into an 
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increase in erosion in 2003 of around 35 % compared to the previous 
period (Table 6). 


In the period from 2004 to 2018, TES and SJN showed lower 
volumes of erosion, possibly related to the typical rapid recovery of 
tropical vegetation in the upper and middle areas (Figure 4), but the 
agricultural area also expands in the form of patches along the two sub- 


basins (Figure 2). 


Cycles of substitution of primary vegetation by areas of agricultural 
and/or livestock use are observed, which may be temporarily abandoned 
to natural regeneration and later reopened for grazing or agriculture 
(Alavez-Vargas et al., 2021). 


In this study, it has been possible to relate the changes in sediment 
production with the cover only when the vegetation cover is modified in 
extensive areas of low slope, such as the increase of grasslands in the low 
areas or when it happens in the steep areas of basin headwaters, such as 
deforestation for agricultural use. A study in the La Antigua basin 
effectively documented a greater production of sediments in grasslands 


compared to forest areas (Martínez et a/., 2009). 


In these case studies, it was found that the topography of the 
Papaloapan is dominated by areas with slopes lower than 15 %, so that, 
while the headwaters can reach steep slopes of more than 2,000 masl, 
more than half of the landscape of each sub-basin consists of wide alluvial 
plains that work as sediment deposit areas (Verstraeten et a/., 2009; Wohl 


et al., 2012). In other words, regardless of their source, the sediments 
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produced are not exported outside the sub-basin, since no statistically 
significant trends were observed in the volume of sediments transported 
by the river between 1986 and 1998 (Alavez-Vargas et al., 2021). 
Similarly, the meta-analysis carried out by Guzha et al. (2018) in the 
African tropics revealed low correlations between forest cover and runoff, 
thus emphasizing the need to improve long-term basin monitoring 


programs to increase the understanding of hydrological responses. 


Potential erosion and historical data scarcity 


The annual production of sediments estimated with MUSLE uses runoff as 
the energy factor of particle detachment and transport, instead of the 
erosivity factor of the precipitation of USLE (Williams € Bernat, 1977). For 
this reason, the availability of observed data for calibration is essential, 
since flow adjustment automatically implies an adjustment to the 
estimated erosion, and its omission can lead to significant under- or over- 


estimates. 


The lack of historical records of erosion measurements was foreseen 
as a limitation to the scope of this study, since it would limit the validation 


of the estimates. Therefore, we proceeded as follows: for the last 
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simulated period (2012-2018), we used the available soil erosion 
cartography of Mexico (INEGI, 2014), which was developed using 
different techniques for the analysis of satellite images, as well as 
information on soil profiles, which shows four categories of erosion based 
on the observed soil attributes (Bolaños-González et a/., 2016). Under the 
criteria of such information, most of the eroded areas in TES were 
classified as light and moderate, presenting a greater extension in the 
middle and lower zones, where their zones of "strong" erosion (bedrock 
exposure in 50-90 % of the erosion polygon) coincide with our high 
erosion results (50-200 ton ha"! year!). 


In SJN, moderate erosion polygons are concentrated in the high and 
low zones in the same areas that we classified as with very high erosion 
(> 200 ton ha! year1?). This comparison is limited by the absence of 
volumetric units in the aforementioned cartography; however, the 
cartography of Borrelli et a/. (2017), estimated at a spatial resolution of 
250 m, identifies numerous points of water erosion greater than 20 ton 
ha"! year “1 in our study area. In contrast, the study conducted by Montes- 
León et al. (2011) classifies most of the Papaloapan basin (and therefore, 
the two sub-basins studied here) in erosion categories above 50 ton ha"! 


year?, 


The divergence between the estimates of the cited studies, together 
with the high level of erosion diagnosed in the country, referred to at the 
beginning of this document, ratifies the importance of having field 


information in order to make a better measurement not only of the current 
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state of soils and vegetation cover of the basins, but also of the scope of 
the results of studies such as this one. Similarly, Tan et al. (2019) 
conclude from their meta-analysis in Southeast Asia that it is urgent to 
aim the efforts towards the in-field generation of good quality data, which 


are essential for modeling and prediction purposes. 


The sources of uncertainty and the scope of modeling 


The SWATCUP software was developed as a calibration tool for SWAT with 
the aim of systematizing, through five possible algorithms, the sensitivity 
and uncertainty analyzes and the parameterization involved in the 
calibration process, given the large number of parameters included in 
each SWAT modeling (Abbaspour, 2015). 


In a comparison of SWATCUP calibration methods conducted on the 
Ganges River, Shivhare, Dikshit and Dwivedi (2018) documented a better 
performance of SUFI-2 over GLUE or ParaSol, supporting its widespread 
use (Abbaspour et al., 2015; Arnold et al., 2012; Chotpantarat 4 
Boonkaewwan, 2018). In this research, the selection of the optimal group 
of parameters for calibration was made based on the: initial 


recommendations of the calibration algorithm and the sensitivity analysis, 
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as well as the guidance of previous studies (Abbaspour et al., 2015; 
Arnold et al., 2012). The optimal set of parameters for each modeling 
included: 1 for surface runoff, 3 for groundwater, 1 for soil water, and 1- 
3 for potential and actual evapotranspiration (Table 3, Table 4, and Table 
5), which have been used for hydrological calibration in other recent 
studies both at a regional scale in Thailand (Chotpantarat 8 
Boonkaewwan, 2018) and at a continental scale in Europe (Abbaspour et 
al., 2015) with good results. Tan et a/. (2019) found good levels of SWAT 
performance in simulating flows in tropical regions from the analysis of 


126 studies from Southeast Asia. 


According to Moriasi et al. (Moriasi et a/., 2007; Moriasi, Gitau, Pai, 
8: Daggupati, 2015) the efficiency values of NSE>0.50 and PBIAS+15 %, 
can be considered satisfactory in the performance of monthly flow 
modelings. Under these criteria, the 12 study periods had fairly good 
percentages of bias (between -8.3 and 7.5), while according to the Nash- 
Sutcliffe efficiency criterion, 6 of 10 calibrated periods were satisfactory, 
the exception being the 1999-2003 period in the two sub-basins, as 
shown in Figure 4d and Figure 5d. The value of the Kling-Gupta objective 
function being close to 1 is just as important as the error of the models 
tending to O, ¡.e., that the highest percentage of the observed data should 
be located in the narrowest prediction interval (95PPU) (R-factor—>.0) of 


the best parameterization (P-factor>1). 


With this combination of criteria, the best simulations of this study 


simultaneously reached KGE and P-factor values greater than 0.70 and 
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correspond to the calibrated periods of 1994-1998, 2004-2011, and 
2012-2018 in TES (Figure 4c, 4d, and 4e) and 1987-1993, 1994-1998, 
and 2012-2018 in SIN (Figure 5b, 5c, and 5e). The performance of the 
1986-1993 simulation in TES is observed in a lower range, with P-factor 
= 0.63 (Figure 4b), while the remaining simulations show errors greater 
than 40 %. As documented by Tan et al. (2019) studies carried out in 
tropical areas under conditions of data scarcity frequently omit the 
description of their error metrics, and these are important to carry out a 
more rigorous evaluation of the results, even if they show, as in this 
study, important anomalies that should be watched for in subsequent 
studies. 


Abbaspour et al. (Abbaspour et a/l., 2007; Abbaspour et al., 2015) 
mention four sources of uncertainty in modeling: 1) conceptual 
simplification, 2) processes not simulated by the model, 3) processes 
simulated by the model, but unknown by the modeler or not modelable 
due to lack of data, and 4) the quality of input data. Points 3 and 4 were 
especially influential in this research. 


Regarding the estimation of water erosion, it was decided to make 
the a priori assumption of the absence of conservation practices (factor 
PUSLE=1 of MUSLE) due to the lack of historical records. Regarding the 
quality of the data, the drastic change in vegetation of the 2003 
scenario—compared to the previous one (1998) and the subsequent one 
(2011)—and the poor performance of the hydrological simulations of 


1999-2003 in the two sub-basins suggest the need to reanalyze that 
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scenario with a different set of images. Possibly, the irregularity of the 
intervals between the evaluated coverage and land use scenarios was a 


source of uncertainty for the model. 


The capacity of the models to reconstruct trajectories of 
hydrological change has as its main limitation the existence (or not) of 
observed historical data to validate the generated scenarios; however, 
these continue to be useful tools to establish baselines (Arheimer € 
Lindstróm, 2019) and, in these cases, the comparison with estimates from 
other models can be used as a strategy for validating the results (Borrelli 
et al., 2017) and limiting the scopes. 


Conclusions 


The response of potential erosion in two medium-scale tropical 
hydrological sub-basins with respect to changes in land use was studied. 
To perform this quantification, the semi-distributed hydrological model 
Soil and Water Assessment Tool (SWAT) was used. 


Of the 12 simulated periods, 2 were invalidated due to anomalies in 
the land use scenario. Based on the combined evaluation of three 
efficiency and two uncertainty statistic data sets, six of the ten remaining 


calibrated periods showed good fitting in the hydrological simulation with 
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KGE>0.70 and P-factor>0.70. According to the estimates made, a 
reduction in the volumes of water erosion was observed in the medium 


term. 


The quality of the results was evaluated through the combined use 
of three efficiency statistics and two measurements of uncertainty. When 
using hydrological models for research and support for public policy 
decision-making, professionals and academicians are recommended to 
select a combination of statistics equivalent to the one used in this study, 
in order to verify jointly and with greater rigor the reliability of the 


modelings. 


The comparison of the volumes of potential erosion estimated with 
respect to the results obtained by other studies in the same study area, 
albeit with different methods and scales, allows validating the results and 


limiting the scopes. 


Likewise, the heterogeneity observed between the studies that were 
reviewed and the high level of erosion diagnosed in the country allows us 
to argue in favor of the need to generate standards and criteria for the 
production of good-quality data for the research and management of 


erosion in Mexico. 


In countries like Mexico, with limited monitoring, the estimation of 
the variables of interest from physical, climate and human data, as well 
as the availability of resampling databases, has popularized the use of 
SWAT. The features of this tool are valuable for prospecting ungauged 


basins or chronosequences of interest; however, its applications for 
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research or management purposes must be supported by error metrics as 
well as by efficiency statistics that clearly measure the scope, limitations 


and anomalies of the results. 


Physically based models are useful tools in tropical regions with a 
scarce historical data to estimate the behavior of hydrological and 
edaphological dynamics, but given the levels of uncertainty found in this 
study, it is considered that they cannot permanently substitute the 


functions of a basin-scale erosion monitoring network. 
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